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ABSTRACT 


Three-finger toxins (TFTs) are well-recognized non- 
enzymatic venom proteins found in snakes. However, 
although TFIs exhibit accelerated evolution, the 
drivers of this evolution remain poorly understood. 
The structural complexes between long-chain 
a-neurotoxins, a subfamily of TFTs, and their nicotinic 
acetylcholine receptor targets have been determined 
in previous research, providing an opportunity to 
address such questions. In the current study, we 
observed several previously identified positively 
selected sites (PSSs) and the highly variable 
C-terminal loop of these toxins at the toxin/receptor 
interface. Of interest, analysis of the molecular 
adaptation of the toxin-recognition regions in the 
corresponding receptors provided no statistical 
evidence for positive selection. However, these 
regions accumulated abundant amino acid variations 
in the receptors from the prey of snakes, suggesting 
that accelerated substitution of TFTs could be a 
consequence of adaptation to these variations. To the 
best of our knowledge, this atypical evolution, initially 
discovered in scorpions, is reported in snake toxins for 
the first time and may be applicable for the evolution 
of toxins from other venomous animals. 
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INTRODUCTION 


Three-finger toxins (TFTs) are among the most abundantly 
secreted and effective components of snake venom. They are 
encoded by a large multigene family and show a diversity of 
functional activities (Fry et al., 2003; Kini, 2011). Members in 
this family contain 57—82 residues and four conserved disulfide 
bridges (Endo & Tamiya, 1987), and are identified by three 
loops extending from a central core resembling their namesake 
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three fingers (Fry et al., 2003). 

The long-term evolutionary process that has resulted in the 
high diversity of TFTs conforms to the birth-and-death model 
of multigene family evolution (Nei et al., 1997). This model of 
evolution produces a series of toxins, allowing snakes to adapt 
to a variety of prey and predators (Nei et al., 1997). Snakes 
employ their venom as a weapon to disable prey (primary 
function) and as a defensive tool against predators (secondary 
function) (Kang et al., 2011). Snake toxins and prey are 
likely involved in a co-evolutionary arms race, whereby evolving 
toxin resistance in prey species and novel toxin evolution in 
snake species exert mutual selective effects (Arbuckle et al., 
2017; Calvete, 2017; Jansa & Voss, 2011; Voss & Jansa, 
2012). However, despite evidence of accelerated evolution 
in the TFT family (Sunagar et al., 2013), the cause of this 
evolution remains vague. 

The TFT family members can target various receptors and 
ion channels with high affinity and specificity (Doley et al., 
2009; Kini, 2011). The o-neurotoxins (a-NTXs) of TFTs 
can interact with nicotinic acetylcholine receptors (nAChRs), 
inhibiting postsynaptic membrane ion flow and leading to 
flaccid paralysis (Barber et al., 2013; Chiappinelli et al., 1996). 
Based on the chain length and number of disulphide bridges, 
a-NTXs are usually divided into short or long a-NTXs. Short 
a-NTXs contain 60-62 amino acid residues with four disulfide 
bonds, whereas long a-NTXs contain 66-75 residues and five 
disulfide bonds (Dajas-Bailador et al., 1998). 

nAChRs are pentameric transmembrane proteins of 
ligand-gated ion channels and are formed by different 
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combinations of five subunits, that is, a, B, y, 6 and e (Corringer 
et al., 2000; Karlin, 2002). Each subunit is composed of a large 
N-terminal extracellular domain, which serves as the major 
binding site for the toxins, followed by four transmembrane 
helices and a small C-terminal extracellular domain (Dellisanti 
et al., 2007; Karlin, 2002). nAChRs can be further divided into 
muscular or neuronal types (Corringer et al., 2000; Wang et 
al., 2002). a-NTXs can potently antagonize the a1 subunit of 
heteropentameric muscle nAChRs ((a@1)28176 in fetal muscle 
and (a1)oB1e6 in adult muscle) and the a7, a8, a9 or a10 
subunits of homopentameric neuronal nAChRs (Karlin, 2002; 
Sine et al., 2013). Crystal structures in the complexes between 
long a-NTXs and related receptors have been identified and 
offer opportunities to explore the molecular mechanism driving 
the accelerated evolution of these toxins (Bourne et al., 2005; 
Dellisanti et al., 2007; Huang et al., 2013). 

In this work, we found several sites previously identified as 
positively selected sites (PSSs) as well as the highly variable 
C-terminal loop of the toxins to be located at the toxin/receptor 


binding interface (Bourne et al., 2005; Sunagar et al., 2013). 


Structural and evolutionary analyses of the toxin-recognition 
region from the receptors (a1, a7, a9 and «10 subunits from 
nAChRs) uncovered an atypical evolution between snakes and 
their prey, in which the amino acid diversity of the nAChR 
toxin-binding regions appeared to drive the adaptive evolution 
of the TFT family. This study showed good agreement with 
our previous research on scorpion toxins and sodium channels 
from their competitors (Zhang et al., 2015). Furthermore, this 
paper provides a broader vision into the evolution between 
venomous animals and their prey/predators. 


MATERIALS AND METHODS 


Sequence analysis 

ClustalX (http://www.clustal.org) was used to align all 
nucleotide and amino acid sequences. A total of 130 long 
a-NTX sequences were aligned and used for sequence logo 
analysis with the WebLogo program (Supplementary Figure 
S1). For the receptors, 76 sequences of muscle-type a1 
nAChRs (three from Reptilia, three from Amphibian, four from 
Aves, 31 from small mammals, and 35 from fishes), 59 
neuronal-tyope a7 nAChRs (one from Gastropoda, two from 
Arthropoda, three from Reptilia, three from Amphibian, three 
from Aves, 18 from small mammals, and 29 from fishes), 
68 neuronal-type a9 nAChRs (three from Reptilia, three from 
Amphibian, four from Aves, 30 from small mammals, and 28 
from fishes), and 52 neuronal-typoe a10 nAChRs (two from 
Reptilia, two from Amphibian, four from Aves, 24 from small 
mammals, and 20 from fishes) (Supplementary Figures S2—S5) 
were also aligned. For positive selection analyses, aligned 
nucleotide sequences of the receptors were used to construct 
neighbor-joining trees. 


WebLogo and ConSurf analyses 


WebLogo can be used to generate sequence logos, which 
are graphical representations of the patterns within multiple 
sequence alignments. The alignments of the aforementioned 
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toxins were used to create sequence logos to identify the 
conservation of each position (Supplementary Figure $1). 
Each logo comprises stacks of letters, one stack for each 
position in the sequence (Crooks et al., 2004). The overall 
height of each stack shows the sequence conservation of that 
position, whereas the height of the symbols within the stack 
indicates the comparative frequency of the amino acid at the 
position (Crooks et al., 2004). Generally, sequence logos 
provide a richer and more precise description of conserved 
and variable regions within sequences. We used the ConSurf 
program (http://consurf.tau.ac.il/) under default parameters to 
calculate conservation scores of the amino acid sequences of 
the related receptors. ConSurf not only depends on sequence 
alignments but also on phylogenetic trees to identify conserved 
and variable regions (Landau et al., 2005). A Bayesian tree 
was constructed using the corresponding alignments with the 
JTT evolutionary substitution model. One of the advantages of 
ConSurf compared with other methods is that the computation 
of the evolutionary rate is more precise when employing the 
empirical Bayesian or maximum-likelinood methods. 


Positive selection analysis 

Excess nonsynonymous substitutions compared with 
synonymous substitutions (w=d\/ds>1) is an important 
sign of positive selection at the molecular level (Yang, 1998; 
Zhu & Gao, 2016). To perform such analysis, we compared 
two pairs of site models (Mla (neutral)/M2a (selection) and 
M7 (beta)/M8 (beta & w)) to measure the selective pressure 
of the receptors to which the long a-NTXs bind (a1, a7, a9 
and «10 subunits from nAChRs). Model M2a and M8 add a 
site class to Mia and M7, respectively, with the free œ ratio 
calculated from the data and used to determine the probability 
of positive selection (Anisimova et al., 2001; Anisimova et al., 
2003; Wong et al., 2004; Yang & Swanson, 2002). As Mia 
and M7 are nested within their respective alternative models 
(M2a and M8) and have two more parameters, x? distribution 
can be used for the likelihood ratio test to compare the fit 
of the two competing models. We used the Bayes Empirical 
Bayes method to calculate the posterior possibility that each 
codon is from the site class of positive selection. The Bayes 
Empirical Bayes method is an improvement of the previous 
Naive Empirical Bayes method and accounts for sampling 
errors in the maximum-likelihood estimates of parameters in 
the model (Nielsen & Yang, 1998). Sites with a high possibility 
(>95%) of coming from the class with œ>1 are likely under 
positive selection and can be analyzed further (Yang, 1998). 


RESULTS 


PSSs of toxins locating on the toxin-receptor complex 
interface 

The N-terminal extracellular domains of nAChRs, which consist 
of a 10 stranded B-sandwich and an N-terminal a-helix, act 
as the major binding sites for long a-NTXs (Changeux et al., 
1970). The three regions of long a-NTXs comprising fingers | 
and Il and the C-terminus are involved in interactions with the 
receptors, with finger Il being the main stabilizing interaction 


center (Bourne et al., 2005). The o211 structure (mouse 
nAChR a1 subunit (PDB entry 2qc1)) can be seen as a 
representative of the N-terminal extracellular domain of the 
nAChR subunit (Dellisanti et al., 2007), in which loops B4-B5 
(loop A), B7-B8 (loop B) and B9-B10 (loop C) serve as the 
principal ligand-binding interfaces (Brejc et al., 2001; Unwin, 
2005) and loop C is the most important region for high affinity 
with long a-NTXs, as revealed by site-directed mutations 
(Fruchart-Gaillard et al., 2002; Levandoski et al., 1999). Loop C 
of «211 is enveloped by fingers | and II and the C-terminal loop 
of the toxin, with finger Il inserted into the ligand-binding site 
wrapped by loops A, B and C of @211 (Figure 1A) (Dellisanti 
et al., 2007). In addition, finger | is sandwiched by loop C, 
whereas finger Ill weakly contributes to «211 binding (Dellisanti 
et al., 2007). 
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Figure 1 «-bungarotoxin interacts with «211 and the sequence 
logos of long a-NTXs 

A: Interaction between a-bungarotoxin (long a-NTXs) and œ211. 
The PSSs 


of the long a-NTXs are in pink. The three loops of #211, major toxin-binding 


a-bungarotoxin is shown in wheat and «211 in light blue. 
regions, are highlighted in blue. B: Sequence logos of long a-NTXs. The 
C-terminal loop of the long a-NTXs involved in receptor binding is indicated 


in a red box. 


We mapped the PSSs of long a-NTXs identified by Sunagar 


et al. (2013) on the complex structure of a-bungarotoxin and 
its receptor (Figure 1A). The toxin-receptor complex model 
showed that some PSSs of the long a-NTXs were located 
on the toxin binding surface with the receptors (Bourne et al., 
2005; Dellisanti et al., 2007; Huang et al., 2013). Almost 
all PSSs in finger Il (Ala31, Phe32, Ser34, Ser35 and 
Val39) of a-bungarotoxin are involved in the interactions with 
the receptors (Dellisanti et al., 2007; Dimitropoulos et al., 
2011; Huang et al., 2013). The same situation appears in 
a-cobratoxin (long a-NTX) since their fast evolved sites (Ala28, 
Phe29, Ser31, lle32 and Arg36) in finger Il also overlap with 
the toxin binding sites of a-cobratoxins (Bourne et al., 2005; 
Dimitropoulos et al., 2011). 

The C-terminal loop of the long a-NTXs also contributes 
to the interaction with related receptors. Multiple sequence 
alignments of long a-NTXs, generated by ClustalX, were used 
to create sequence logos (Figure 1B, Supplementary Figure 
S1) (Crooks et al., 2004). The C-terminal loop of the long 
a-NTXs, indicated in the red box in Figure 1B, was highly 
variable based on the WebLogo analyses, and the whole C-tail 
involves extensive insertions and deletions. Thus, this loop 
might undergo positive selection despite the technical difficulty 
in detecting PSSs from indel-containing sequences. 


No evidence for positive selection in the toxin-recognition 
regions of receptors 


Snakes employ their venom to immobilize various prey, 
including snails, insects, fishes, toads, lizards, chickens, small 
mammals and even other snakes (Kang et al., 2011). The 
maximum-likelihnood models of codon substitutions were used 
to identify selective pressure in the toxin-recognition regions 
of the receptors from the prey of snakes. However, no 
positive selection signals were detected in the receptors of long 
a-NTXs (a1, «7, a9 and a10 subunits) (Table 1, Supplementary 
Tables S1—S4). The maximum-likelinood estimates under MO 
showed that the average œ ratios for all receptor sequence 
pairs ranged from 0.02 to 0.07. M2a and M8 detected no 
evolution-accelerating sites. The ws under M8 of the a1 and 
&œ10 subunits of nAChRs from snake prey equaled 1 (Table 
1). Under the M8 model, the «7 and a9 subunits of nAChRs 
showed ws>1 (Table 1), but their proportions (p1) equaled O, 
proving that no PSSs existed (Supplementary Tables 81-84). 


Table 1 Parameter estimates and likelihood ratio statistics (2^1) for different subunit types of nAChRs 


Type l (M1a) l (M2a) We (M2a) 
a1 subunit —12874.0 -—12874.0 1.00 
a7 subunit —13163.9 -—13163.9 1.00 
a9 subunit —11856.3 -11856.3 1.00 
a10 subunit —9954.2 —9954.2 1.00 


2Al 1 (M7) 1(M8) ws (M8) 2A 


O -12617.4 —-12607.4 1.00 20 
O -12785.7 -12785.7 34.66 0 
O -11626.9 -11626.9 1.42 0 
0 -9746.8 -—9746.8 1.00 0 


l is the log likelihood; 2A/ is between null models and their alternative models: M1a/M2a and M7/M8. Two o values in M2a and M8 are 


boldfaced. 
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Snake toxins bind to variable regions of their receptors 
Although PSSs were detected in the toxins, none were found in 
the toxin-recognition regions of the involved receptors. Based 
on the complex models between the long a-NTXs and related 
receptors, we further analyzed the evolutionary conservation 
of the N-terminal extracellular domain regions of the nAChR «a1 
and a7 subunits from snake prey using ConSurf (Dellisanti et 
al., 2007; Huang et al., 2013) (Figure 2). Our results indicated 
that loop C demonstrated the greatest variation among the 
three receptor loops (Figure 2). 
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Figure 2 ConSurf results of x-bungarotoxin with the receptors 
ConSurf results of a-bungarotoxin with «1 (PDB entry 2qc1) and a7 (PDB 
entry 4hqp) nAChRs are shown in A and B, respectively. ConSurf provides 
evolutionary conservation scores for nAChRs. PSSs of finger | and finger II 
and the residues of the C-terminal loop in a-bungarotoxin are highlighted in 
red (PSSs of finger II are partly hidden). Toxin-binding regions (loop A, loop 
B and loop C) of the a1 and a7 nAChRs are spherized. Binding sites variable 


in loop C are indicated. 
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We analyzed the interaction pairs in the a-bungarotoxin 
and a1 nAChR subunit complex, which were lle11-Phe189/ 
Proi94, Val31-Tyr93/Asp99/Phe100, Phe32-Tyr93/Phe100/ 
Trp149/Tyr190, Val39-Val188/Tyr190, His68-Pro194, Pro69- 
Ser191, Lys70-Pro194 and Gln71-Pro194, with the toxin sites 
corresponding to the PSSs in fingers | and Il and the C-terminal 
loop (Dellisanti et al., 2007; Dimitropoulos et al., 2011). The 
loop A (Tyr93, Asp99 and Phe100), loop B (Trp149) and loop 
C (Val188, Phe189, Tyri90, Ser191, Pro194 and Pro197) 
were involved in the interaction with the PSSs of the toxin 
binding surface. Although Tyr93 in loop A and Trp149 in loop B 
were conserved, Asp99 and Phe100 in loop A exhibited some 
variability (99: Asp/Asn and 100: Phe/Tyr). In contrast, the 
binding sites in loop C were more variable (188: Val/Arg/Lys, 
189: Phe/Tyr/Thr, 191: Ala/Ser/Thr/Pro and 194: Pro/Gln) 
(Figure 2A). 

For the a-bungarotoxin and a7 nAChR subunit complex, 
the interaction pairs include lle11-Phe183/Lys188, Phe32- 
Tyr91/Trpo145/Arg182/Tyr184, Val39-Arg182/Tyr184, His68- 
Lys188, Pro69-Glu185, Lys70-Cys186/Lys188 and Gln71- 
Lys188, with the sites of the toxins also corresponding to 
the PSSs in fingers | and Il and the C-terminal loop (Huang 
et al., 2013). Compared with the conserved residues in loop 
A and loop B (Tyr91 and Trp145), the binding sites in loop 
C were diverse (182: Lys/Arg/Leu/Ser/Asn, 183: Phe/Tyr/lle, 
185: Glu/Asp/Asn/Gly and 188: Lys/Asp/Glu/Pro/Gln) (Figure 
2B). The additional PSSs (Val31, Ser34 and Ser35) in finger 
Il contacted residues of the complementary subunit, and 
included Ser32, Ser34, Leu36, Tro53, Glnd55, Gln114, Leu116 
and Asp160 (Huang et al., 2013). These residues in the 
complementary subunit may provide the driving force for the 
additional sites in finger Il. 

Taken together, our results suggest that the variable 
toxin-recognition region in the receptors might drive the 
accelerated evolution of the toxin-binding residues. 


DISCUSSION 


By mapping the PSSs of the long a-NTXs on the toxin-receptor 
complex, we found that several of them are located on the 
toxin-binding surface of the receptor (Figure 1A). In addition 
to these PSSs, high sequence diversity was also observed 
in the C-terminal loop of the long a-NTXs. Thus, given its 
importance in the interaction with receptors, we surmised that 
it could be an accelerated substituted region for adaptation to 
receptor variability (Figure 1B). Several PSSs were detected 
in finger IIl of the long a-NTXs, although finger III contributed 
little to the binding. This may be due to its high flexibility 
in the complex. Other mechanisms may be involved in the 
accelerated substitutions of this finger, which requires further 
investigation. 

We further observed the amino acid variability of the 
principal toxin-recognition regions (mainly in loop C) in related 
nAChR subunits (Figure 2). Compared with loop A and loop B, 
loop C exhibited the greatest variation due to its predominant 
role in the toxin-receptor interactions. Thus, we concluded that 
the evolutionary variability of the toxin-recognition regions of 


the receptors is a possible driver for the accelerated evolution 
of the toxins. 

Short-chain a-NTXs can bind to nAChRs with high affinity 
(Tremeau et al., 1995). Loop II of short-chain a-NTXs pushes 
into the ligand-binding pocket of nAChRs, whereas the tips 
of loops | and Ill contact nAChRs only in a ‘surface-touch’ 
way (Mordvintsev et al., 2005). Previous study on Torpedo 
californica showed that the C-loop is vital for the binding of 
short a-NTXs to nAChRs (Bourne et al., 2005; Mordvintsev et 
al., 2005). However, as there are no specific coordinate files of 
complexes between short a-NTXs and their targets, hindering 
further study. 

Our previous study on the evolution of scorpion toxins and 
voltage-gated sodium (Nay) channels indicated that variability 
of the toxin-recognition regions in the Nay channels from 
scorpion predators and prey is a putative driver of the 
accelerated evolution of the functional regions of scorpion 
toxins following gene duplication (Zhang et al., 2015). Similarly, 
PSSs exist in the binding interface of the long a-NTXs, 
though only amino acid variability was detected in the 
principal toxin-recognition regions of the related nAChR 
subunits, suggesting that amino-acid substitutions on the 
toxin-recognition surface in related nAChR subunits could 


provide a force driving the accelerated evolution of the toxins. 


Thus, the atypical co-evolutionary manner between snake 
toxins and their receptors is similar to our previous research 
on scorpion toxins and their targets (Zhang et al., 2015). We 
supposed that accelerated evolution of the receptor-bound 
regions of the snake toxins is a consequence of adaptation 
to variable receptors of their prey. From the viewpoint of 
a co-evolutionary arms race between predators and prey 
(Arbuckle et al., 2017; Calvete, 2017; Jansa & Voss, 2011; 
Voss & Jansa, 2012), it appears that prey might exert greater 
selective pressure on their predators, as described in our 
current study. As this evolutionary manner has been shown 
to occur in two distant species, we believe it will be revealed in 
more toxins from diverse venomous animals. 
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